#!/usr/bin/env python2.6
# -*- coding: utf-8 -*-
# Copyright (c) 2008 Qtrac Ltd. All rights reserved.
# This program or module is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version. It is provided for educational
# purposes and is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.


import sys
import pickle
import string
import csv
import os
import io
import numpy as np
from scipy import stats
import time


"""Trimming the dataraw file to be used for the entire procedure"""


def get_epi_dataraw(data_file_raw,path,partition,traits,genelist):
#independent script    
    """This script trims the tab-delimited file to matrix form"""
    """All monomorph snps are removed"""
    print("36) EPIC Call get_epi_dataraw""")

    epirawfile = path+data_file_raw#job["data_file2"]    
    data_raw=list(open(epirawfile,"r"))#,encoding="utf-8")#only in py3.1

    """Organize and clean-up raw data"""
    """Remove tabs, split (\t)"""

    for subj in range(0,len(data_raw)):
        data_raw[subj]=data_raw[subj].split("\t")
    """..'"' and empty spaces are removed, and in that order"""
    for subj in range(0,len(data_raw)):
        for col in range(0,len(data_raw[0])):
            data_raw[subj][col] = data_raw[subj][col].strip('"')
    for subj in range(0,len(data_raw)):
        for col in range(0,len(data_raw[0])):
            data_raw[subj][col] = data_raw[subj][col].strip()

#Empty strings allocated 0, which is OK for genes, but for traits?
    """Empty cells are filled with 0 for removal of invalid date entries"""
    for subj in range(0,len(data_raw)):
        for col in range(0,len(data_raw[0])):
            if data_raw[subj][col] == '' or data_raw[subj][col] == '"':
               data_raw[subj][col] = '0'
               
    """Create matrix"""

    data_raw_matrix = np.matrix(data_raw)

    """Remove partitions, traits, and snps not included in the analysis, but is present in the raw-file"""
    #Reduce matrix to include selected parttition, traits and genes only
    gene_list_limit = data_raw_matrix.shape[1]-genelist[0]#+2+len(traits)
    gene_in_use = genelist[0:gene_list_limit]   
    index_list = [0]+[partition]+traits+gene_in_use
    data_raw_matrix_trim=data_raw_matrix[0:,index_list]
    partition_use = data_raw_matrix_trim[0,1]
    traits_use = range(2,2+len(traits))
    
    """Exclude monomorphic snps"""
    genelist_first = 2+len(traits)
    genelist_last = len(index_list)
    snp_use = {}#dict of snps tested

    for i in range(genelist_first,genelist_last):
        snp_use[i] = data_raw_matrix_trim[0,i]    
#	print(i)
#	print(snp_use[i])
    gt1 = 0
    gt2 = 0
    gt3 = 0
    snp_summary =[]
    snp_momomorph_all = []
    mono_discard = []#list of itmes (snps) to be deleted from matrix
    monomorf_general = open("Epicalc/MonomorphGen.txt", "w")   
    monomorf_general_header = 'SNP_number' + ' ' + "rs-number"+"\n"
    monomorf_general.write(monomorf_general_header)    

    #number of snps
    numb_snps_tot = genelist_last - genelist_first#only valid if range is used?
    numb_snps_sub = 0#to be used for calculating correct significantlevel
    for snps in range(genelist_first,genelist_last):
        gt1 = 0
        gt2 = 0
        gt3 = 0
        snp_check = data_raw_matrix_trim[1:,[snps]].astype(float)       
        for gt in range(0, len(snp_check)):
            if snp_check[[gt]]==1:
               gt1+=1
            if snp_check[[gt]]==2:
               gt2+=1
            if snp_check[[gt]]==3:
               gt3+=1    
        item_list = list((gt1, gt2, gt3))
        
        gtzero = 0
        mono_list = []
        for snp in range(0,len(item_list)):                     
            if item_list[snp] == 0:
               gtzero +=1
        if gtzero == 2:
           mono_list = [snps,snp_use.get(snps)]
           snp_momomorph_all.append(mono_list)
           mono_discard.append(snps)

           line = "{0} {1}\n".format(mono_list[0],mono_list[1]) 
           monomorf_general.write(line)
           numb_snps_sub+=1
    #Number of snps included
    numb_snps = numb_snps_tot - numb_snps_sub
    snps_in_file = [numb_snps_tot,numb_snps,numb_snps_sub]

    """Monomophic snps are removed from matrix"""
    data_raw_matrix_trimmed = np.delete(data_raw_matrix_trim, mono_discard, 1)
    print("End data trim")
    return data_raw_matrix_trimmed,snps_in_file,partition_use,traits_use


def epi_selection(job):
    print("37) Executing calculations")

    data_raw_matrix_trimmed = job['data_file_trimmed']  
#    print(data_raw_matrix_trimmed)	  
    class_select = float(job['class'][0])

    """Selection of entries with the given class number"""
    class_check = data_raw_matrix_trimmed[1:,1].astype(float)       
    clasno = 0
    subj_select = [0]#zero to include subject id
    for clch in range(0,len(class_check)):
        if class_check[[clch]]==class_select:
           subj_select.append(clch+1)           
           clasno+=1

    """Class-specific matrix"""
    data_raw_matrix_trimmed_class = data_raw_matrix_trimmed.take(subj_select,axis=0)

    subj_numb = data_raw_matrix_trimmed_class.shape[0]-1 
    partition = job['partition_use']#only used for summary
    traitlist = job['traits_use']
    gene_first_trimmed = 2 + len(traitlist)
    gene_last_trimmed = data_raw_matrix_trimmed_class.shape[1]-1

    data_file_raw = job['data_file_raw']#only to print name of data file in Summary
    snps_in_file=job['snps_remaining']

         
    lntrait = len(traitlist)
#    nom_sign_level = 0.05#from config!!!!!
    num_snp_tests = snps_in_file[1]*(snps_in_file[1]-1)/2

#Bonferoni corection
    nom_sign_level = job["significance_level"]
    Bonf_correct = nom_sign_level/(num_snp_tests*lntrait)
    #should number partitions be included???done before
    
    snplist = range(gene_first_trimmed,gene_last_trimmed+1)
    
    """Det vil særdeles hensigtsmæssigt at køre databas-struktur!!!"""

    """Summary file"""
    

    summary_job = open("Summary", "w")
    summary_output1 = 'Start of calculations'+' '+str(time.asctime())+'\n\n'
    summary_job.write(summary_output1)
    summary_output0 = "Summary of entire currated data file"+"\n\n"

    summary_job.write(summary_output0)

    summary_output2 = "File                       "+ str(data_file_raw.split("/")[-1]) +"\n\n"
    summary_job.write(summary_output2)

    summary_output9 = "Partition                 "+ str(partition) +"\n"
    summary_job.write(summary_output9)

    summary_output10 = "Class                     "+ str(class_select) +"\n"
    summary_job.write(summary_output10)

    summary_output13 = "Number of subjects        "+ str(subj_numb) +"\n\n"
    summary_job.write(summary_output13)


    summary_output4 = "Total number of SNPs       "+ str(snps_in_file[0]) +"\n"
    summary_job.write(summary_output4)  
    summary_output5 = "Number of monomorph SNPs   "+ str(snps_in_file[2]) +"\n"
    summary_job.write(summary_output5)  
    summary_output6 = "Number of valid SNPs       "+ str(snps_in_file[1]) +"\n\n"
    summary_job.write(summary_output6)  

    summary_output16 = "Maximum number of possible interactions   "+ str(num_snp_tests) +"\n\n"
    summary_job.write(summary_output16)  


    summary_output14 = "Nominal signficance level "+ str(nom_sign_level) +"\n"
    summary_job.write(summary_output14)  

    summary_output15 = "Bonferoni correction      "+ str(Bonf_correct) +"\n\n"
    summary_job.write(summary_output15)  

    if clasno == 0:
       summary_output15 = "No subjects in the class"+"\n"
       summary_job.write(summary_output15)  
       return 


    """Here the for-loops are defined for traits if more than one and for all genes"""  

    print("class "+str(class_select))
    print("Trait name")
    for tr in range(0,len(traitlist)):
        trait = traitlist[tr]
        name_trait = data_raw_matrix_trimmed_class[0,traitlist[tr]]


        print(name_trait)
        trait_use = []
        epi_sum_out = []
        """Calculating mean and variance for the entire file"""
        """This should be moved as it is common to all"""
        for subj in range(1,len(data_raw_matrix_trimmed)):
            subj_val = float(data_raw_matrix_trimmed[subj,trait])  
            if subj_val > 0:#zero filled in in missing entries
               trait_use.append(subj_val)
               
        trait_mean = np.average(trait_use)
        trait_var = np.var(trait_use)
        trait_mean_class =  np.average(data_raw_matrix_trimmed_class[1:,traitlist[tr]].astype(float))
#        print(trait_mean_class)   
        trait_var_class =  np.var(data_raw_matrix_trimmed_class[1:,traitlist[tr]].astype(float))
#        print(trait_var_class)   
      
        
        
        #Allel frequence file
        genotype_allefr = open("GenotypesAllelfr "+name_trait, "w")
        genotype_allefr_header = 'SNP1' +' '+'AA'+ ' ' + 'Aa'+ ' '+ 'aa' + ' ' + "Allel_A" + ' ' 'Allel_a'\
                              +' '+'SNP2' +' '+'BB'+ ' ' + 'Bb'+ ' '+ 'bb' + ' ' + "Allel_B" + ' ' 'Allel_b'+ "\n"
        genotype_allefr.write(genotype_allefr_header)

        #Effect file
        genotype_effects = open("GenotypeEffects "+name_trait, "w")
        genotype_effects_header = 'SNP1' +' '+'Add_A'+ ' ' + 'Add_a'+ ' '+ 'DomAA'\
          + ' ' +'DomAa' + ' ' + 'Domaa' + ' ' 'SNP2' + ' ' + "Add_B" + ' ' 'Add_b'\
          +' '+'DomBB' +' '+'DomBb'+ ' ' + 'Dombb'+' '+ 'Interactions'+' ' \
          + 'AddAdd_AB'+' '+'AddAdd_Ab'+' '+'AddAdd_aB'+' '+'AddAdd_ab'+ ' '\
          +'AddDom_AxBB'+' '+'AddDom_AxBb' + ' ' + 'AddDom_Axbb' + ' '\
          +'AddDom_axBB'+' '+'AddDom_axBb'+' '+'AddDom_axbb'+' '\
          +'AddDom_AABx'+' '+'AddDom_AaBx' + ' ' + 'AddDom_aaBx' + ' '\
          +'AddDom_AAbx'+' '+'AddDom_Aabx'+' '+'AddDom_aabx'+ ' '\
          +'DomDomAABB'+' '+'DomDomAABb'+' '+'DomDomAAbb'+' '+'DomDomAaBB'+' '+'DomDomAaBa'+' '\
          +'DomDomAabb'+' '+'DomDomaaBB'+' '+'DomDomaaBb'+' '+'DomDomaabb'+"\n"
        genotype_effects.write(genotype_effects_header)    

        #Variance detail file
        genotype_variance = open("GenotypeVariance "+name_trait, "w")
        genotype_variance_header = 'SNP1' +' '+'AddVarAa'+ ' ' + 'DomVarAA'+ ' '+ 'DomVarAa'\
          + ' ' +'DomVaraa' + ' ' + 'SNP2' + ' ' + "AddVarBb" + ' ' 'DomVarBB'\
          +' '+'DomVarBb' +' '+'DomVarbb'+ ' ' + 'Interactions'+' '+'AddAddVarAB'+' ' \
          + 'AddAddVarAb'+' '+'AddAddVaraB'+' '+'AddAddVarab'+ ' '\
          +'AddDomVar_AxBB'+' '+'AddDomVar_AxBb' + ' ' + 'AddDomVar_Axbb' + ' '\
          +'AddDomVar_axBB'+' '+'AddDomVar_axBb'+' '+'AddDomVar_axbb'+' '\
          +'AddDomVar_AABx'+' '+'AddDomVar_AaBx' + ' ' + 'AddDomVar_aaBx' + ' '\
          +'AddDomVar_AAbx'+' '+'AddDomVar_Aabx'+' '+'AddDomVar_aabx'+ ' '\
          +'DomDomVarAABB'+' '+'DomDomVarAABb'+' '+'DomDomVarAAbb'+' '+'DomDomVarAaBB'+' '\
          +'DomDomVarAaBa'+' '+'DomDomVarAabb'+' '+'DomDomVaraaBB'+' '+'DomDomVaraaBb'+' '\
          +'DomDomaabb'+"\n"
        genotype_variance.write(genotype_variance_header)    
        #Variance summary file
        genotype_variance_sum = open("GenotypeVarSum"+name_trait, "w")
        genotype_variance_header_sum = 'SNP1' +' '+'SNP2'+ ' ' + 'PhenoVar'+ ' '\
            +'GenoVar'+' '+'Additive'+' '+'Dominant'+' '+'AddAdd'+' '+'AddDom(A)'+' '\
            +'AddDom(B)'+' '+'DomDom'+' '+'EpiVar'+' '+'RelVar'+' '+'Geno/Pheno'\
            +' '+'Epi/Pheno'+' '+'Epi/Geno'+' '+'Check'+' '+'F-stat'\
            + ' '+ 'dfdenom'+' '+'p-values'+"\n"
        genotype_variance_sum.write(genotype_variance_header_sum)  

        """bonf is defined above, this is just for test purposes"""  

        for snpA in range(0,len(snplist)):
            snp1 = snplist[snpA]
            name_snpA = data_raw_matrix_trimmed_class[0,snplist[snpA]]
            for snpB in range((snpA+1),len(snplist)):
                snp2 = snplist[snpB]
                name_snpB = data_raw_matrix_trimmed_class[0,snplist[snpB]]
                
                """Created matrix with selected trait and snps"""
                data_use_matrix = get_selection(data_raw_matrix_trimmed_class,trait,snp1,snp2)
                name_SNP1 = data_use_matrix[0,1]
                name_SNP2 = data_use_matrix[0,2]

                """Calculate allel frequences"""
                allel_val= get_allelfreqences(data_use_matrix,trait,snp1,snp2,name_SNP1,name_SNP2,Bonf_correct)
#return case_valid, allelfreq, gtall
                """Check for HWE; skip if non-equilibirum"""
                if allel_val[0] == 0:
                   continue
                #monomorph is not registered, but can be if momomorph
                #is run for the partion separately. Do!!, saves time!! Use above modified
                #should be done as the first call to griepei when jobs are created. 
                """Calculating means"""
                geno_mean = epistasis_calc_mean(data_use_matrix,allel_val[1])
#return ltraitgt, mtraitgt, vtraitgt, TotalgMeans
                
                """Calculating significans of epistasis"""
                episign = epivar_sign(geno_mean[0],geno_mean[1],geno_mean[2])
#return Fstat, dfdenom, Fpvalue
                """Check for epistasis significance; skip if non-significant"""
                if episign[2] >= Bonf_correct:
                    continue
                
                """Calculating total variance"""
                tot_var = epistasis_calc_totvar(allel_val[1],geno_mean[1],geno_mean[3],name_SNP1,name_SNP2)
                                     
#return PhenVar, TGenVar
                
                """Calculating epistasis"""
                
                epi_calc = epistasis_calc_var(allel_val[1],geno_mean[1],geno_mean[3],episign,tot_var,name_SNP1,name_SNP2)

#return gteffects, gtvariance, gtvariance_sum

                """Output"""
                genotype_allefr.write(allel_val[2])
                genotype_effects.write(epi_calc[0])
                genotype_variance.write(epi_calc[1])
                genotype_variance_sum.write(epi_calc[2])
                epi_sum_out.append(epi_calc[2][0])#no header!!

#        print(epi_sum_out)
	epistasis_significant = len(epi_sum_out)
        relative_sign = float(epistasis_significant)/float(num_snp_tests)

        summary_output3 = "Trait: "+ name_trait +"\n\n"
        summary_job.write(summary_output3)  
        summary_output7 = "  Phenotypic mean            "+ str(trait_mean) +"\n"
        summary_job.write(summary_output7)
        summary_output8 = "  Phenotypic variance        "+ str(trait_var) +"\n\n"
        summary_job.write(summary_output8)
        
        summary_output11 = "  Phenotypic mean class "+str(class_select)+"\t"   + str(trait_mean_class) +"\n"
        summary_job.write(summary_output11)

        summary_output12 = "  Phenotypic variance class "+ str(trait_var_class) +"\n\n"
        summary_job.write(summary_output12)
        
        summary_output17 = "  Significant interactions "+ str(epistasis_significant)+" ("+ str(relative_sign)+")" +"\n\n"
        summary_job.write(summary_output17)

        
def get_selection(filename,trait,snp1,snp2):
 #   print("get_selection")
    data_use_matrix = np.take(filename, (trait,snp1,snp2),1)
    """Deletes subjects if one or both genetypes are missing or trait value is missing"""
    subj_discard = []
    for subj in range(1,len(filename)):
        subj_var = [float(filename[subj,trait]),float(filename[subj,snp1]),float(filename[subj,snp2])]  
        if np.prod(subj_var) == 0:
           subj_discard.append(subj)

    data_use_matrix = np.delete(data_use_matrix, subj_discard, 0)

    return data_use_matrix

def get_allelfreqences(data_use_matrix,trait,snp1,snp2,name_SNP1,name_SNP2,Bonf_correct):
#    print("get_allelfreqences")
    sum_cases = data_use_matrix.shape[0]-1
    gt1 = 0
    gt2 = 0
    gt3 = 0
    gt_countA = data_use_matrix[1:,[1]].astype(int)
    for gt in range(0, len(gt_countA)):
        if gt_countA[[gt]]==1:
           gt1+=1
        if gt_countA[[gt]]==2:
           gt2+=1
        if gt_countA[[gt]]==3:
           gt3+=1    
    genotypeA = list((gt1, gt2, gt3))
    
    gt1 = 0
    gt2 = 0
    gt3 = 0
    gt_countB = data_use_matrix[1:,[2]].astype(int)
    for gt in range(0, len(gt_countB)):
        if gt_countB[[gt]]==1:
           gt1+=1
        if gt_countB[[gt]]==2:
           gt2+=1
        if gt_countB[[gt]]==3:
           gt3+=1    
    genotypeB = list((gt1, gt2, gt3))

    """Check for monomorph snps. If so, further calculation is skipped"""
     
    gtzeroA = 0
    gtzeroB = 0 
    for snpA in range(0,len(genotypeA)):                     
        if genotypeA[snpA] == 0:
           gtzeroA +=1
    if gtzeroA == 2:
       case_valid = 0
       return case_valid, name_SNP1# comma alone ok but may be confusing
    for snpB in range(0,len(genotypeB)):                     
        if genotypeB[snpB] == 0:
           gtzeroB +=1
    if gtzeroB == 2:
       case_valid = 0
       return case_valid, name_SNP2
    
    """Allele frequencies are calculated"""
    allelA = float(2*genotypeA[0] + genotypeA[1])/(sum_cases*2)
    allela = float(2*genotypeA[2] + genotypeA[1])/(sum_cases*2)
    allelB = float(2*genotypeB[0] + genotypeB[1])/(sum_cases*2)
    allelb = float(2*genotypeB[2] + genotypeB[1])/(sum_cases*2)
    allelfreq = list((allelA, allela, allelB, allelb))

    """HWE calculations"""
    #Gene A
    obs_sum_geneA = sum(genotypeA)
    exp_geneA = [allelfreq[0]*allelfreq[0]*obs_sum_geneA,
                2*allelfreq[0]*allelfreq[1]*obs_sum_geneA,
                allelfreq[1]*allelfreq[1]*obs_sum_geneA]

    obs_exp_diff_powA = [pow(genotypeA[0]-exp_geneA[0],2),
                         pow(genotypeA[1]-exp_geneA[1],2),
                         pow(genotypeA[2]-exp_geneA[2],2)]
    chis_valA = [obs_exp_diff_powA[0]/exp_geneA[0],
                 obs_exp_diff_powA[1]/exp_geneA[1],
                 obs_exp_diff_powA[2]/exp_geneA[2]]
    chis_val_sumA = sum(chis_valA)
    hwe_geneA = stats.chisqprob(chis_val_sumA,2)

    #Gene B
    obs_sum_geneB = sum(genotypeB)
    exp_geneB = [allelfreq[2]*allelfreq[2]*obs_sum_geneB,
                2*allelfreq[2]*allelfreq[3]*obs_sum_geneB,
                allelfreq[3]*allelfreq[3]*obs_sum_geneB]

    obs_exp_diff_powB = [pow(genotypeB[0]-exp_geneB[0],2),
                         pow(genotypeB[1]-exp_geneB[1],2),
                         pow(genotypeB[2]-exp_geneB[2],2)]
    chis_valB = [obs_exp_diff_powB[0]/exp_geneB[0],
                 obs_exp_diff_powB[1]/exp_geneB[1],
                 obs_exp_diff_powB[2]/exp_geneB[2]]
    chis_val_sumB = sum(chis_valB)
    hwe_geneB = stats.chisqprob(chis_val_sumB,2)

    gtall = ""
    if hwe_geneA >= Bonf_correct and hwe_geneB >= Bonf_correct:
        gtall = "{0} {1} {2} {3} {4} {5} {6} {7} {8} {9} {10} {11}\n".format(name_SNP1,genotypeA[0],
           genotypeA[1],genotypeA[2],
           allelA,allela,name_SNP2,genotypeB[0],genotypeB[1],genotypeB[2],
           allelB,allelb)
        case_valid = 1
    else:
        case_valid = 0
        
    return case_valid, allelfreq, gtall
##FLYT epistasis-sign hertil så vi vi ikke spilder tid med detaljerede bergeninger


def epistasis_calc_mean(data_use_matrix,allelfreq):
 #   print("epistasis_calc_mean")

    #Overall mean
    all_mean = np.average(data_use_matrix[1:,0].astype(float))
    #To be printed to file containing all the means and variances of 2-snps results
#    print(all_mean)
    """Traivalues for every combinations of two snps"""
    traitgt11 = []
    traitgt12 = []
    traitgt13 = []
    traitgt21 = []
    traitgt22 = []
    traitgt23 = []
    traitgt31 = []
    traitgt32 = []
    traitgt33 = []

    """Calculating the trait means for 2-snp interactions"""
    
    for gt in range(1, len(data_use_matrix)):
        if data_use_matrix[gt,1].astype(int)==1 and data_use_matrix[gt,2].astype(int)==1:
           traitgt11.append(data_use_matrix[gt,0].astype(float))
        if data_use_matrix[gt,1].astype(int)==1 and data_use_matrix[gt,2].astype(int)==2:
           traitgt12.append(data_use_matrix[gt,0].astype(float))
        if data_use_matrix[gt,1].astype(int)==1 and data_use_matrix[gt,2].astype(int)==3:
           traitgt13.append(data_use_matrix[gt,0].astype(float))
        if data_use_matrix[gt,1].astype(int)==2 and data_use_matrix[gt,2].astype(int)==1:
           traitgt21.append(data_use_matrix[gt,0].astype(float))
        if data_use_matrix[gt,1].astype(int)==2 and data_use_matrix[gt,2].astype(int)==2:
           traitgt22.append(data_use_matrix[gt,0].astype(float))
        if data_use_matrix[gt,1].astype(int)==2 and data_use_matrix[gt,2].astype(int)==3:
           traitgt23.append(data_use_matrix[gt,0].astype(float))
        if data_use_matrix[gt,1].astype(int)==3 and data_use_matrix[gt,2].astype(int)==1:
           traitgt31.append(data_use_matrix[gt,0].astype(float))
        if data_use_matrix[gt,1].astype(int)==3 and data_use_matrix[gt,2].astype(int)==2:
           traitgt32.append(data_use_matrix[gt,0].astype(float))
        if data_use_matrix[gt,1].astype(int)==3 and data_use_matrix[gt,2].astype(int)==3:
           traitgt33.append(data_use_matrix[gt,0].astype(float))

    """Number of subjects in each group"""
#    global ltraitgt
    ltraitgt = [len(traitgt11),len(traitgt12),len(traitgt13),len(traitgt21),\
                len(traitgt22),len(traitgt23),len(traitgt31),len(traitgt32),\
                len(traitgt33)]
    """Average of traits, trait-menas"""
    mtraitgt11 = np.average(traitgt11)
    mtraitgt12 = np.average(traitgt12)
    mtraitgt13 = np.average(traitgt13)
    mtraitgt21 = np.average(traitgt21)
    mtraitgt22 = np.average(traitgt22)
    mtraitgt23 = np.average(traitgt23)
    mtraitgt31 = np.average(traitgt31)
    mtraitgt32 = np.average(traitgt32)
    mtraitgt33 = np.average(traitgt33)

    mtraitgt = [mtraitgt11,mtraitgt12,mtraitgt13,mtraitgt21,\
                mtraitgt22,mtraitgt23,mtraitgt31,mtraitgt32,mtraitgt33]
    #NB traits without content (0) is nan
    for i in range(len(ltraitgt)):
        if ltraitgt[i] == 0 :mtraitgt[i]=0
#    global vtraitgt
    vtraitgt11 = np.var(traitgt11)
    vtraitgt12 = np.var(traitgt12)
    vtraitgt13 = np.var(traitgt13)
    vtraitgt21 = np.var(traitgt21)
    vtraitgt22 = np.var(traitgt22)
    vtraitgt23 = np.var(traitgt23)
    vtraitgt31 = np.var(traitgt31)
    vtraitgt32 = np.var(traitgt32)
    vtraitgt33 = np.var(traitgt33)
    vtraitgt = [vtraitgt11,vtraitgt12,vtraitgt13,vtraitgt21,\
                vtraitgt22,vtraitgt23,vtraitgt31,vtraitgt32,vtraitgt33]

    """Mean values conditioned on HWE of snps"""

  #Mean values for genotypes calculated as HWE frequences = OK
    #gMean is eq. to gVaacc etc in R-script
    gMean11 = allelfreq[0]*allelfreq[0]*allelfreq[2]*allelfreq[2]*mtraitgt[0]
    gMean12 = 2*allelfreq[0]*allelfreq[0]*allelfreq[2]*allelfreq[3]*mtraitgt[1]
    gMean13 = allelfreq[0]*allelfreq[0]*allelfreq[3]*allelfreq[3]*mtraitgt[2]
    gMean21 = 2*allelfreq[0]*allelfreq[1]*allelfreq[2]*allelfreq[2]*mtraitgt[3]
    gMean22 = 4*allelfreq[0]*allelfreq[1]*allelfreq[2]*allelfreq[3]*mtraitgt[4]
    gMean23 = 2*allelfreq[0]*allelfreq[1]*allelfreq[3]*allelfreq[3]*mtraitgt[5]
    gMean31 = allelfreq[1]*allelfreq[1]*allelfreq[2]*allelfreq[2]*mtraitgt[6]
    gMean32 = 2*allelfreq[1]*allelfreq[1]*allelfreq[2]*allelfreq[3]*mtraitgt[7]
    gMean33 = allelfreq[1]*allelfreq[1]*allelfreq[3]*allelfreq[3]*mtraitgt[8]

    gMean = [gMean11,gMean12,gMean13,gMean21,gMean22,gMean23,gMean31,gMean32,gMean33]
    #TotalgMeans = TotgValues in original R-script
    TotalgMeans = np.sum(gMean)

    return ltraitgt, mtraitgt, vtraitgt, TotalgMeans

def epistasis_calc_totvar(allelfreq,mtraitgt,TotalgMeans,name_SNP1,name_SNP2):

    #NB change to AABB etc 
    gVar11 = allelfreq[0]*allelfreq[0]*allelfreq[2]*allelfreq[2]*pow(mtraitgt[0],2)
    gVar12 = 2*allelfreq[0]*allelfreq[0]*allelfreq[2]*allelfreq[3]*pow(mtraitgt[1],2)
    gVar13 = allelfreq[0]*allelfreq[0]*allelfreq[3]*allelfreq[3]*pow(mtraitgt[2],2)
    gVar21 = 2*allelfreq[0]*allelfreq[1]*allelfreq[2]*allelfreq[2]*pow(mtraitgt[3],2)
    gVar22 = 4*allelfreq[0]*allelfreq[1]*allelfreq[2]*allelfreq[3]*pow(mtraitgt[4],2)
    gVar23 = 2*allelfreq[0]*allelfreq[1]*allelfreq[3]*allelfreq[3]*pow(mtraitgt[5],2)
    gVar31 = allelfreq[1]*allelfreq[1]*allelfreq[2]*allelfreq[2]*pow(mtraitgt[6],2)
    gVar32 = 2*allelfreq[1]*allelfreq[1]*allelfreq[2]*allelfreq[3]*pow(mtraitgt[7],2)
    gVar33 = allelfreq[1]*allelfreq[1]*allelfreq[3]*allelfreq[3]*pow(mtraitgt[8],2)

    PhenVar = np.sum([gVar11,gVar12,gVar13,gVar21,gVar22,gVar23,gVar31,gVar32,gVar33])
    #For checksum: TGenVar the same as in R-script
    TGenVar = PhenVar - pow(TotalgMeans,2)

    return PhenVar, TGenVar
    
def epistasis_calc_var(allelfreq,mtraitgt,TotalgMeans,episign,tot_var,name_SNP1,name_SNP2):


    """Genotypic values"""

#Allele A OK, fix A
    addValAABB = allelfreq[0]*allelfreq[2]*allelfreq[2]*mtraitgt[0]
    addValAABb = 2*allelfreq[0]*allelfreq[2]*allelfreq[3]*mtraitgt[1]
    addValAAbb = allelfreq[0]*allelfreq[3]*allelfreq[3]*mtraitgt[2]
    addValAaBB = allelfreq[1]*allelfreq[2]*allelfreq[2]*mtraitgt[3]
    addValAaBb = 2*allelfreq[1]*allelfreq[2]*allelfreq[3]*mtraitgt[4]
    addValAabb = allelfreq[1]*allelfreq[3]*allelfreq[3]*mtraitgt[5]

    addValA = [addValAABB,addValAABb,addValAAbb,addValAaBB,addValAaBb,addValAabb]
    addValATot = np.sum(addValA)
 
#Allel a OK, fix a 
    addValAaBB = allelfreq[0]*allelfreq[2]*allelfreq[2]*mtraitgt[3]
    addValAaBb = 2*allelfreq[0]*allelfreq[2]*allelfreq[3]*mtraitgt[4]
    addValAabb = allelfreq[0]*allelfreq[3]*allelfreq[3]*mtraitgt[5]
    addValaaBB = allelfreq[1]*allelfreq[2]*allelfreq[2]*mtraitgt[6]
    addValaaBb = 2*allelfreq[1]*allelfreq[2]*allelfreq[3]*mtraitgt[7]
    addValaabb = allelfreq[1]*allelfreq[3]*allelfreq[3]*mtraitgt[8]

    addVala = [addValAaBB,addValAaBb,addValAabb,addValaaBB,addValaaBb,addValaabb]
    addValaTot = np.sum(addVala)

#Allel B OK, fix B
    addValAABB = allelfreq[0]*allelfreq[0]*allelfreq[2]*mtraitgt[0]
    addValAaBB = 2*allelfreq[0]*allelfreq[1]*allelfreq[2]*mtraitgt[3]
    addValaaBB = allelfreq[1]*allelfreq[1]*allelfreq[2]*mtraitgt[6]
    addValAABb = allelfreq[0]*allelfreq[0]*allelfreq[3]*mtraitgt[1]
    addValAaBb = 2*allelfreq[0]*allelfreq[1]*allelfreq[3]*mtraitgt[4]
    addValaaBb = allelfreq[1]*allelfreq[1]*allelfreq[3]*mtraitgt[7]

    addValB = [addValAABB,addValAaBB,addValaaBB,addValAABb,addValAaBb,addValaaBb]
    addValBTot = np.sum(addValB)
 
#Allel b OK, fix b
    addValAABb = allelfreq[0]*allelfreq[0]*allelfreq[2]*mtraitgt[1]
    addValAaBb = 2*allelfreq[0]*allelfreq[1]*allelfreq[2]*mtraitgt[4]
    addValaaBb = allelfreq[1]*allelfreq[1]*allelfreq[2]*mtraitgt[7]
    addValAAbb = allelfreq[0]*allelfreq[0]*allelfreq[3]*mtraitgt[2]
    addValAabb = 2*allelfreq[0]*allelfreq[1]*allelfreq[3]*mtraitgt[5]
    addValaabb = allelfreq[1]*allelfreq[1]*allelfreq[3]*mtraitgt[8]

    addValb = [addValAABb,addValAaBb,addValaaBb,addValAAbb,addValAabb,addValaabb]
    addValbTot = np.sum(addValb)
  
    """Additive effects, zero centered OK"""

    addEffA = addValATot - TotalgMeans
    addEffa = addValaTot - TotalgMeans
    addEffB = addValBTot - TotalgMeans
    addEffb = addValbTot - TotalgMeans
    #Could this be used to say something about increasing/decreasing the trait?
    addCheckSNP1 = [addEffA,addEffa]
    addCheckSNP2 = [addEffB,addEffb]

    """Additive variances"""
    AddVarSNP1 = 2*(allelfreq[0]*pow(addEffA,2) + allelfreq[1]*pow(addEffa,2))
    AddVarSNP2 = 2*(allelfreq[2]*pow(addEffB,2) + allelfreq[3]*pow(addEffb,2))

    """Total additive variance OK"""
    AddVarTot = AddVarSNP1 + AddVarSNP2

    """Dominant values, effects and variance"""

#Conditonal genotypic values SNP1 OK
    domValAA = pow(allelfreq[2],2)*mtraitgt[0] + 2*allelfreq[2]*allelfreq[3]*mtraitgt[1] + pow(allelfreq[3],2)*mtraitgt[2]
    domValAa = pow(allelfreq[2],2)*mtraitgt[3] + 2*allelfreq[2]*allelfreq[3]*mtraitgt[4] + pow(allelfreq[3],2)*mtraitgt[5]
    domValaa = pow(allelfreq[2],2)*mtraitgt[6] + 2*allelfreq[2]*allelfreq[3]*mtraitgt[7] + pow(allelfreq[3],2)*mtraitgt[8]
#Conditonal genotypic values SNP2 OK
    domValBB = pow(allelfreq[0],2)*mtraitgt[0] + 2*allelfreq[0]*allelfreq[1]*mtraitgt[3] + pow(allelfreq[1],2)*mtraitgt[6]
    domValBb = pow(allelfreq[0],2)*mtraitgt[1] + 2*allelfreq[0]*allelfreq[1]*mtraitgt[4] + pow(allelfreq[1],2)*mtraitgt[7]
    domValbb = pow(allelfreq[0],2)*mtraitgt[2] + 2*allelfreq[0]*allelfreq[1]*mtraitgt[5] + pow(allelfreq[1],2)*mtraitgt[8]

#Dominant effects OK
    domEffAA = domValAA - TotalgMeans - addEffA - addEffA
    domEffAa = domValAa - TotalgMeans - addEffA - addEffa
    domEffaa = domValaa - TotalgMeans - addEffa - addEffa
    domEffBB = domValBB - TotalgMeans - addEffB - addEffB
    domEffBb = domValBb - TotalgMeans - addEffB - addEffb
    domEffbb = domValbb - TotalgMeans - addEffb - addEffb
    domCheckSNP1 = [domEffAA,domEffAa,domEffaa]
    domCheckSNP2 = [domEffBB,domEffBb,domEffbb]
    
#Dominant variance weigthed OK
    domVarAA = allelfreq[0]*allelfreq[0]*pow(domEffAA,2)
    domVarAa = 2*allelfreq[0]*allelfreq[1]*pow(domEffAa,2)
    domVaraa = allelfreq[1]*allelfreq[1]*pow(domEffaa,2)
    domVarBB = allelfreq[2]*allelfreq[2]*pow(domEffBB,2)
    domVarBb = 2*allelfreq[2]*allelfreq[3]*pow(domEffBb,2)
    domVarbb = allelfreq[3]*allelfreq[3]*pow(domEffbb,2)
    vardomCheckSNP1 = [domVarAA,domVarAa,domVaraa]
    vardomCheckSNP2 = [domVarBB,domVarBb,domVarbb]

    domVarTot = np.sum(vardomCheckSNP1 + vardomCheckSNP2)

    """Additive*additive values, effects and variance"""

#Additive*Additive values OK

    addaddValAB = allelfreq[0]*allelfreq[2]*mtraitgt[0]+allelfreq[0]*allelfreq[3]*mtraitgt[1]\
                  +allelfreq[1]*allelfreq[2]*mtraitgt[3]+allelfreq[1]*allelfreq[3]*mtraitgt[4]
    addaddValAb = allelfreq[0]*allelfreq[2]*mtraitgt[1]+allelfreq[0]*allelfreq[3]*mtraitgt[2]\
                  +allelfreq[1]*allelfreq[2]*mtraitgt[4]+allelfreq[1]*allelfreq[3]*mtraitgt[5]
    addaddValaB = allelfreq[0]*allelfreq[2]*mtraitgt[3]+allelfreq[0]*allelfreq[3]*mtraitgt[4]\
                  +allelfreq[1]*allelfreq[2]*mtraitgt[6]+allelfreq[1]*allelfreq[3]*mtraitgt[7]
    addaddValab = allelfreq[0]*allelfreq[2]*mtraitgt[4]+allelfreq[0]*allelfreq[3]*mtraitgt[5]\
                  +allelfreq[1]*allelfreq[2]*mtraitgt[7]+allelfreq[1]*allelfreq[3]*mtraitgt[8]

    addaddCheckGene = [addaddValAB,addaddValAb,addaddValaB,addaddValab]

#Additive*Additive effects OK
    addaddEffAB = addaddValAB - TotalgMeans - addEffA - addEffB
    addaddEffAb = addaddValAb - TotalgMeans - addEffA - addEffb
    addaddEffaB = addaddValaB - TotalgMeans - addEffa - addEffB
    addaddEffab = addaddValab - TotalgMeans - addEffa - addEffb
    
    addaddEff = [addaddEffAB,addaddEffAb,addaddEffaB,addaddEffab]
 
#Additive*Additive variance OK

    addaddVarAB = allelfreq[0]*allelfreq[2]*pow(addaddEffAB,2)
    addaddVarAb = allelfreq[0]*allelfreq[3]*pow(addaddEffAb,2)
    addaddVaraB = allelfreq[1]*allelfreq[2]*pow(addaddEffaB,2)
    addaddVarab = allelfreq[1]*allelfreq[3]*pow(addaddEffab,2)

    addaddVar = [addaddVarAB,addaddVarAb,addaddVaraB,addaddVarab]
    addaddVarTot = 4*np.sum(addaddVar)

    """Additive*dominance values, effects and variances"""

#Additive*dominance values OK
    #SNP1
    adddomValAxBB = allelfreq[0]*mtraitgt[0] + allelfreq[1]*mtraitgt[3]
    adddomValAxBb = allelfreq[0]*mtraitgt[1] + allelfreq[1]*mtraitgt[4]
    adddomValAxbb = allelfreq[0]*mtraitgt[2] + allelfreq[1]*mtraitgt[5]
    adddomValaxBB = allelfreq[0]*mtraitgt[3] + allelfreq[1]*mtraitgt[6]
    adddomValaxBb = allelfreq[0]*mtraitgt[4] + allelfreq[1]*mtraitgt[7]
    adddomValaxbb = allelfreq[0]*mtraitgt[5] + allelfreq[1]*mtraitgt[8]

    adddomValA = [adddomValAxBB,adddomValAxBb,adddomValAxbb,\
                  adddomValaxBB,adddomValaxBb,adddomValaxbb]
 
    #SNP2
    adddomValAABx = allelfreq[2]*mtraitgt[0] + allelfreq[3]*mtraitgt[1]
    adddomValAaBx = allelfreq[2]*mtraitgt[3] + allelfreq[3]*mtraitgt[4]    
    adddomValaaBx = allelfreq[2]*mtraitgt[6] + allelfreq[3]*mtraitgt[7]
    adddomValAAbx = allelfreq[2]*mtraitgt[1] + allelfreq[3]*mtraitgt[2]
    adddomValAabx = allelfreq[2]*mtraitgt[4] + allelfreq[3]*mtraitgt[5]
    adddomValaabx = allelfreq[2]*mtraitgt[7] + allelfreq[3]*mtraitgt[8]

    adddomValB = [adddomValAABx,adddomValAaBx,adddomValaaBx,\
                  adddomValAAbx,adddomValAabx,adddomValaabx]
    
#Additive*dominance effects OK
    #SNP1
    adddomEffAxBB = adddomValAxBB - TotalgMeans - addEffA - addEffB - addEffB\
                    -domEffBB - addaddEffAB - addaddEffAB
    adddomEffAxBb = adddomValAxBb - TotalgMeans - addEffA - addEffB - addEffb\
                    -domEffBb - addaddEffAB - addaddEffAb
    adddomEffAxbb = adddomValAxbb - TotalgMeans - addEffA - addEffb - addEffb\
                    -domEffbb - addaddEffAb - addaddEffAb
    adddomEffaxBB = adddomValaxBB - TotalgMeans - addEffa - addEffB - addEffB\
                    -domEffBB - addaddEffaB - addaddEffaB
    adddomEffaxBb = adddomValaxBb - TotalgMeans - addEffa - addEffB - addEffb\
                    -domEffBb - addaddEffaB - addaddEffab
    adddomEffaxbb = adddomValaxbb - TotalgMeans - addEffa - addEffb - addEffb\
                    -domEffbb - addaddEffab - addaddEffab

    adddomEffA = [adddomEffAxBB,adddomEffAxBb,adddomEffAxbb,adddomEffaxBB,adddomEffaxBb,adddomEffaxbb]
 
    #SNP2
    adddomEffAABx = adddomValAABx - TotalgMeans - addEffA - addEffA - addEffB\
                    -domEffAA - addaddEffAB - addaddEffAB
    adddomEffAaBx = adddomValAaBx - TotalgMeans - addEffA - addEffa - addEffB\
                    -domEffAa - addaddEffAB - addaddEffaB
    adddomEffaaBx = adddomValaaBx - TotalgMeans - addEffa - addEffa - addEffB\
                    -domEffaa - addaddEffaB - addaddEffaB
    adddomEffAAbx = adddomValAAbx - TotalgMeans - addEffA - addEffA - addEffb\
                    -domEffAA - addaddEffAb - addaddEffAb
    adddomEffAabx = adddomValAabx - TotalgMeans - addEffA - addEffa - addEffb\
                    -domEffAa - addaddEffAb - addaddEffab
    adddomEffaabx = adddomValaabx - TotalgMeans - addEffa - addEffa - addEffb\
                    -domEffaa - addaddEffab - addaddEffab

    adddomEffB = [adddomEffAABx,adddomEffAaBx,adddomEffaaBx,adddomEffAAbx,adddomEffAabx,adddomEffaabx]

#Additive*dominance variance OK
    #SNP1
    addomVarAxBB = allelfreq[0]*allelfreq[2]*allelfreq[2]*pow(adddomEffAxBB,2)
    addomVarAxBb = 2*allelfreq[0]*allelfreq[2]*allelfreq[3]*pow(adddomEffAxBb,2)
    addomVarAxbb = allelfreq[0]*allelfreq[3]*allelfreq[3]*pow(adddomEffAxbb,2)
    addomVaraxBB = allelfreq[1]*allelfreq[2]*allelfreq[2]*pow(adddomEffaxBB,2)
    addomVaraxBb = 2*allelfreq[1]*allelfreq[2]*allelfreq[3]*pow(adddomEffaxBb,2)
    addomVaraxbb = allelfreq[1]*allelfreq[3]*allelfreq[3]*pow(adddomEffaxbb,2)

    adddomVarA = [addomVarAxBB,addomVarAxBb,addomVarAxbb,addomVaraxBB,addomVaraxBb,addomVaraxbb]
 
    #SNP2
    addomVarAABx = allelfreq[0]*allelfreq[0]*allelfreq[2]*pow(adddomEffAABx,2)
    addomVarAaBx = 2*allelfreq[0]*allelfreq[1]*allelfreq[2]*pow(adddomEffAaBx,2)
    addomVaraaBx = allelfreq[1]*allelfreq[1]*allelfreq[2]*pow(adddomEffaaBx,2)
    addomVarAAbx = allelfreq[0]*allelfreq[0]*allelfreq[3]*pow(adddomEffAAbx,2)
    addomVarAabx = 2*allelfreq[0]*allelfreq[1]*allelfreq[3]*pow(adddomEffAabx,2)
    addomVaraabx = allelfreq[1]*allelfreq[1]*allelfreq[3]*pow(adddomEffaabx,2)
        
    adddomVarB = [addomVarAABx,addomVarAaBx,addomVaraaBx,addomVarAAbx,addomVarAabx,addomVaraabx]

    adddomVarTotA = 2*np.sum(adddomVarA)
    adddomVarTotB = 2*np.sum(adddomVarB)
    adddomVarTot = adddomVarTotA+adddomVarTotB 
    
##Dominant*dominance effects and variances (no intermediary values)

#Dominance*dominance effects OK

    domdomEffAABB = mtraitgt[0]-TotalgMeans-addEffA-addEffA-addEffB-addEffB\
            -domEffAA-domEffBB-addaddEffAB-addaddEffAB-addaddEffAB-addaddEffAB\
            -adddomEffAxBB-adddomEffAxBB-adddomEffAABx-adddomEffAABx
    domdomEffAABb = mtraitgt[1]-TotalgMeans-addEffA-addEffA-addEffB-addEffb\
            -domEffAA-domEffBb-addaddEffAB-addaddEffAB-addaddEffAb-addaddEffAb\
            -adddomEffAxBb-adddomEffAxBb-adddomEffAABx-adddomEffAAbx
    domdomEffAAbb = mtraitgt[2]-TotalgMeans-addEffA-addEffA-addEffb-addEffb\
            -domEffAA-domEffbb-addaddEffAb-addaddEffAb-addaddEffAb-addaddEffAb\
            -adddomEffAxbb-adddomEffAxbb-adddomEffAAbx-adddomEffAAbx
    domdomEffAaBB = mtraitgt[3]-TotalgMeans-addEffA-addEffa-addEffB-addEffB\
            -domEffAa-domEffBB-addaddEffAB-addaddEffAB-addaddEffaB-addaddEffaB\
            -adddomEffAxBB-adddomEffaxBB-adddomEffAaBx-adddomEffAaBx
    domdomEffAaBb = mtraitgt[4]-TotalgMeans-addEffA-addEffa-addEffB-addEffb\
            -domEffAa-domEffBb-addaddEffAB-addaddEffAb-addaddEffaB-addaddEffab\
            -adddomEffAxBb-adddomEffaxBb-adddomEffAaBx-adddomEffAabx
    domdomEffAabb = mtraitgt[5]-TotalgMeans-addEffA-addEffa-addEffb-addEffb\
            -domEffAa-domEffbb-addaddEffAb-addaddEffAb-addaddEffab-addaddEffab\
            -adddomEffAxbb-adddomEffaxbb-adddomEffAabx-adddomEffAabx
    domdomEffaaBB = mtraitgt[6]-TotalgMeans-addEffa-addEffa-addEffB-addEffB\
            -domEffaa-domEffBB-addaddEffaB-addaddEffaB-addaddEffaB-addaddEffaB\
            -adddomEffaxBB-adddomEffaxBB-adddomEffaaBx-adddomEffaaBx
    domdomEffaaBb = mtraitgt[7]-TotalgMeans-addEffa-addEffa-addEffB-addEffb\
            -domEffaa-domEffBb-addaddEffaB-addaddEffaB-addaddEffab-addaddEffab\
            -adddomEffaxBb-adddomEffaxBb-adddomEffaaBx-adddomEffaabx
    domdomEffaabb = mtraitgt[8]-TotalgMeans-addEffa-addEffa-addEffb-addEffb\
            -domEffaa-domEffbb-addaddEffab-addaddEffab-addaddEffab-addaddEffab\
            -adddomEffaxbb-adddomEffaxbb-adddomEffaabx-adddomEffaabx


    domdomEff = [domdomEffAABB,domdomEffAABb,domdomEffAAbb,domdomEffAaBB,\
                 domdomEffAaBb,domdomEffAabb,domdomEffaaBB,domdomEffaaBb,domdomEffaabb]
    
#Dominance*dominance variance OK
    domdomVarAABB = allelfreq[0]*allelfreq[0]*allelfreq[2]*allelfreq[2]*pow(domdomEffAABB,2)
    domdomVarAABb = 2*allelfreq[0]*allelfreq[0]*allelfreq[2]*allelfreq[3]*pow(domdomEffAABb,2)
    domdomVarAAbb = allelfreq[0]*allelfreq[0]*allelfreq[3]*allelfreq[3]*pow(domdomEffAAbb,2)
    domdomVarAaBB = 2*allelfreq[0]*allelfreq[1]*allelfreq[2]*allelfreq[2]*pow(domdomEffAaBB,2)
    domdomVarAaBb = 4*allelfreq[0]*allelfreq[1]*allelfreq[2]*allelfreq[3]*pow(domdomEffAaBb,2)
    domdomVarAabb = 2*allelfreq[0]*allelfreq[1]*allelfreq[3]*allelfreq[3]*pow(domdomEffAabb,2)
    domdomVaraaBB = allelfreq[1]*allelfreq[1]*allelfreq[2]*allelfreq[2]*pow(domdomEffaaBB,2)
    domdomVaraaBb = 2*allelfreq[1]*allelfreq[1]*allelfreq[2]*allelfreq[3]*pow(domdomEffaaBb,2)
    domdomVaraabb = allelfreq[1]*allelfreq[1]*allelfreq[3]*allelfreq[3]*pow(domdomEffaabb,2)

    
    domdomVar = [domdomVarAABB,domdomVarAABb,domdomVarAAbb,domdomVarAaBB,\
                 domdomVarAaBb,domdomVarAabb,domdomVaraaBB,domdomVaraaBb,domdomVaraabb]
    domdomVarTot = np.sum(domdomVar)

    gteffects = "{0} {1} {2} {3} {4} {5} {6} {7} {8} {9} {10} {11}  \
{12} {13} {14} {15} {16} {17} {18} {19} {20} {21} {22} {23} {24} {25} \
{26} {27} {28} {29} {30} {31} {32} {33} {34} {35} {36}\n".format(name_SNP1,addEffA,
           addEffa,domEffAA,domEffAa,domEffaa,name_SNP2,addEffB,addEffb,
           domEffBB,domEffBb,domEffbb,addaddEffAB,addaddEffAb,addaddEffaB,addaddEffab,
           adddomEffA[0],adddomEffA[1],adddomEffA[2],adddomEffA[3],adddomEffA[4],adddomEffA[5],
           adddomEffB[0],adddomEffB[1],adddomEffB[2],adddomEffB[3],adddomEffB[4],adddomEffB[5],
           domdomEff[0],domdomEff[1],domdomEff[2],domdomEff[3],domdomEff[4],
           domdomEff[5],domdomEff[6],domdomEff[7],domdomEff[8])
    
    TotGenVar = np.sum(AddVarTot+domVarTot+addaddVarTot+adddomVarTot+domdomVarTot)

    CheckGenVar = tot_var[1] - TotGenVar#kan droppe når konfidence

    gtvariance = "{0} {1} {2} {3} {4} {5} {6} {7} {8} {9}  {10} {11} {12} {13} \
{14} {15} {16} {17} {18} {19} {20} {21} {22} {23} {24} {25} \
{26} {27} {28} {29} {30} {31} {32} {33} {34}\n".format(name_SNP1,
           AddVarSNP1,domVarAA,domVarAa,domVaraa,name_SNP2,AddVarSNP2,
           domVarBB,domVarBb,domVarbb,addaddVarAB,addaddVarAb,addaddVaraB,
           addaddVarab,addomVarAxBB,addomVarAxBb,addomVarAxbb,addomVaraxBB,
           addomVaraxBb,addomVaraxbb,addomVarAABx,addomVarAaBx,addomVaraaBx,
           addomVarAAbx,addomVarAabx,addomVaraabx,domdomVarAABB,domdomVarAABb,
           domdomVarAAbb,domdomVarAaBB,domdomVarAaBb,domdomVarAabb,domdomVaraaBB,
           domdomVaraaBb,domdomVaraabb)
     
 #Check can be deleted when we are confident!!   

    """Summary of variances"""
    epivar = np.sum(addaddVarTot+adddomVarTotA+adddomVarTotB+domdomVarTot)
    relGP = tot_var[1]/tot_var[0]
    relEP = epivar/tot_var[0]
    relEG = epivar/tot_var[1]
    
    gtvariance_sum = "{0} {1} {2} {3} {4} {5} {6} {7} {8} {9} {10}  {11} {12} {13} \
{14} {15} {16} {17}\n".format(name_SNP1,name_SNP2,tot_var[0],tot_var[1],
           AddVarTot,domVarTot,addaddVarTot,adddomVarTotA,adddomVarTotB,                                   
           domdomVarTot,epivar,relGP,relEP,relEG,CheckGenVar,episign[0],
           episign[1],episign[2])

    return gteffects, gtvariance, gtvariance_sum

def epivar_sign(ltraitgt,mtraitgt,vtraitgt):
#significance tests
    #mtraitgt mean of groups
    """Average of simple genotypic means"""
    gmaverage = np.average(mtraitgt)
    """Calculation of non-epistatic values"""

    gAA = np.average([mtraitgt[0],mtraitgt[1],mtraitgt[2]])
    gAa = np.average([mtraitgt[3],mtraitgt[4],mtraitgt[5]])
    gaa = np.average([mtraitgt[6],mtraitgt[7],mtraitgt[8]])
    gBB = np.average([mtraitgt[0],mtraitgt[3],mtraitgt[6]])
    gBb = np.average([mtraitgt[1],mtraitgt[4],mtraitgt[7]])
    gbb = np.average([mtraitgt[2],mtraitgt[5],mtraitgt[8]])

    """Non-epistatic genotypic values """

    neAABB = gAA + gBB -  gmaverage
    neAABb = gAA + gBb -  gmaverage
    neAAbb = gAA + gbb -  gmaverage
    neAaBB = gAa + gBB -  gmaverage
    neAaBb = gAa + gBb -  gmaverage
    neAabb = gAa + gbb -  gmaverage
    neaaBB = gaa + gBB -  gmaverage
    neaaBb = gaa + gBb -  gmaverage
    neaabb = gaa + gbb -  gmaverage
    nevect = [neAABB,neAABb,neAAbb,neAaBB,neAaBb,neAabb,neaaBB,neaaBb,neaabb]

    """Calculating significance"""

    evect = []
    varmean =[]
    for i in range(0,8):
        eve = ltraitgt[i]*pow(mtraitgt[i]- nevect[i],2)
        evect.append(eve)
        varm = ltraitgt[i]*vtraitgt[i]
        varmean.append(varm)

    EMS = np.sum(evect)/4
    RMS = np.sum(varmean)/sum(ltraitgt)
    dfdenom = sum(ltraitgt)-9
    if dfdenom <= 0:
       Fstat = "Less than 10 cases"
       dfdenom = ""
       Fpvalue = ""
       return Fstat, dfdenom, Fpvalue
    if RMS == 0:#only if var is zero in all groups, which is a rare, but possible event
       Fstat = "F-statistics cannot be calculated"
       dfdenom = ""
       Fpvalue = ""
       return Fstat, dfdenom, Fpvalue

    Fstat = EMS/RMS
    """F-stat significance = stats.fprob(4,dfdenom,Fstat)
       From scipy.special.fdtrc is (df nomo,df denom, F-value)"""
    Fpvalue = stats.fprob(4,dfdenom,Fstat)
    return Fstat, dfdenom, Fpvalue

    
def get_epi_datarawc():
    epirawfile = "data/M101 comp  140210W.csv"
    xxx=csv.reader(open(epirawfile,"rb"))
    print(epirawfile)
    print(xxx)
    type(xxx)
    csvfile = open(epirawfile, "rb").read()

    dialect = csv.Sniffer().sniff(epirawfile)
    print(csvfile)
    print("number of items 1")
#epi_partition()
    
#epirawfile,xxx = get_epi_dataraw()
#get_epi_dataraw()

##this is ok, drop write.
